Purcell's Swimmer

We make a self-propelling Purcell Swimmer and dip it in a highly viscous fluid.

import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
import seaborn as sns; sns.set(); sns.set_style("whitegrid")
plt.rcParams["animation.html"] = "jshtml"
from celluloid import Camera
import stokeslets as slts
import purcell as pur
import matplotlib.patches as ptc

figfile = "figs/"

Purcell's Description

Purcell imagines a swimmer made of three 'links' or rods connected by two joints. The geometry of the swimmer can be described at any given point in time by the two angles $\theta_1$ and $\theta_2$. It undergoes a cycle as shown below from $S_1$ to $S_5$ moving its arms alternately.

Geometry

We need to define the position of the rods and the motors. Let's construct Purcell's and place the end points of the rod in the $S_1$ configuration. The rod is parametrised by $\eta$ — the ratio of the body length to the arm length — and $\phi$ = $\max_t(\theta_1(t)) = \max_t(\theta_2(t))$.

et = 2 #ratio of the body to one arm
ph = np.pi/3
size = 1 #total length

s = np.zeros([4,3]);
a = size/(2+et) #one-arm length
b = a*et #body length
s[1,0:2] = np.array([-b/2,0]) #motor 1
s[2,0:2] = np.array([b/2,0]) #motor 2
s[0,0:2] = np.array([-a*np.cos(ph)-b/2,a*np.sin(ph)]) #arm 1
s[3,0:2] = np.array([a*np.cos(ph)+b/2,a*np.sin(ph)]) #arm 2
r = np.array([1,2]) #motor specifications

fig,ax = plt.subplots(figsize=(4,4))

ax.plot(s.T[0],s.T[1],'k-',alpha=0.5)

arc1 = ptc.Arc((s.T[0,2],s.T[1,2]),0.2,0.2,angle=0,
               theta1=0,theta2=ph*180/np.pi,color='r',linewidth=1.5)

ax.add_patch(arc1)

plt.annotate(r"$\phi$",(s.T[0,2]+0.15,s.T[1,2]+0.075))

plt.xlim(-1,1); plt.ylim(-1,1);
plt.savefig(figfile+"purcell_swimmer_configuration.png",dpi=300,bbox_inches='tight')
plt.show()

Dynamics

Define the time vector.

dt = 0.01
T = 16
tim = np.arange(0,T,dt)

Let's define $\tau_1 = [0,0,q_1(t)]$ and $\tau_2 = [0,0,q_2(t)]$.

tau = np.array([np.zeros(np.shape(tim)+(3,))]*np.size(r))

In Purcell's paper, the swimmer moves one arm at a given time. This can be achieved by a constant torque difference between the two motors. To describe this motion, we need a parameter $P = |q_2(0)|$ which is the maximum torque exerted by one of the motors. For simplicity, we assume that the other motor exerts no torque when the first is switched on.

P = 1.463 #maximum power
tau = pur.stroke_input(P,tim,r,strokes=2) #the second motor exerts no torque when the first is moving

plt.plot(tim,tau[0,:,-1],'-',alpha=0.9);
plt.plot(tim,tau[1,:,-1],'-',alpha=0.9);
plt.legend(["$q_1 (t)$","$q_2 (t)$"],loc="upper right");
plt.xlabel(r"$t \ (a.u)$"); plt.ylabel(r"$\tau_z \ (a.u)$")
plt.savefig(figfile+"torque_over_time.png",dpi=300,bbox_inches="tight")

It Swims!

R = slts.mesher()

_,pswimmer = slts.evolve(tau,tim,R,r,s,k=100,e=0.3,c=0.6)

fig = plt.figure(figsize = (4.5,4.5))
plt.xlim(-1,1); plt.ylim(-1,1)
camera = Camera(fig)
for i in range(0,np.shape(pswimmer)[0],8):
    plt.plot(pswimmer[i].T[0],pswimmer[i].T[1],'k-', linewidth=2, alpha=0.5, label='_nolegend_')
    plt.plot(np.mean(pswimmer[i].T[0]),np.mean(pswimmer[i].T[1]),'ro', linewidth=2, alpha=0.5)
    camera.snap()
plt.legend(["CoM"])
plt.close()

animation = camera.animate()
animation.save(figfile+"purcell_swimmer.gif",writer = 'imagemagick',fps = 20,dpi=300)
animation

Once Loop Reflect

Purcell Strokes

We can try finding the torque difference required for a Purcell stroke by trial and error, but it would be easier to sample different $P$ values and see which corresponds to a Purcell stroke.

Defining $\theta_1$ and $\theta_2$ by Purcell's convention as in the figure above, we can come up with a function that must be equal to zero for a Purcell stroke after half a cycle:

$$\Delta \phi = (\phi-\theta_2)+(\phi+\theta_1)$$

print("P0 = ",pur.delphihalf(tau,pswimmer,ph))

P0 =  0.04077787615948769

Now, we can try sampling over the parameter space to see which values of $P$ give $\Delta\phi_{1/2} = 0$ where $\Delta \phi _{1/2}$ is $\Delta \phi$ after half a cycle is completed.

et = 2 #eta value
ph = np.pi/3 #phi value

res = 0.1
P_vec = np.arange(0.6,1.6+res,res) #range of P values

delphiH = np.zeros(np.shape(P_vec))

for i,P in enumerate(P_vec):
    s,r,tim = pur.init(et,ph,T=9)
    tau = pur.stroke_input(P,tim,r,strokes=9/8)
    R = slts.mesher()
    _,pswimmer = slts.evolve(tau,tim,R,r,s,k=100,e=0.3,c=0.6)
    delphiH[i] = pur.delphihalf(tau,pswimmer,ph)
    print("Progress:",np.round((i+1)/np.size(P_vec)*100,2),"%")

Progress: 9.09 %
Progress: 18.18 %
Progress: 27.27 %
Progress: 36.36 %
Progress: 45.45 %
Progress: 54.55 %
Progress: 63.64 %
Progress: 72.73 %
Progress: 81.82 %
Progress: 90.91 %
Progress: 100.0 %

plt.figure(figsize=(5,3.5))
plt.plot(P_vec,delphiH,'o',color='seagreen');
plt.ylabel(r"$\Delta\phi_{1/2}$")
plt.xlabel(r"$P$")
plt.savefig(figfile+"delphiH_vs_P_example.png",dpi=300,bbox_inches="tight")
plt.show()

The nice one-to-one mapping simplifies our job by quite a bit. We can use bisection method to find $P_0$.